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Chapter  1 


Grid  Generation  for  DNS/LES 


1.1  Two-Dimensional  Grid  Generation 

An  elliptic  grid  generation  method  first  proposed  by  Spekrense  (1995)  is  used  to  generate  2D 
grids.  The  elliptic  grid  generation  method  is  based  on  a  composite  mapping,  which  is  consisted 
of  a  nonlinear  transfinite  algebraic  transformation  and  an  elliptic  transformation.  The  algebraic 
transformation  maps  the  computational  space  C  onto  a  parameter  space  V,  and  the  elliptic  trans¬ 
formation  maps  the  parameter  space  on  to  the  physical  domain  V .  The  computational  space, 
parameter  space,  and  the  physical  domain  are  illustrated  in  Figure  1.1. 
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Figure  1.1:  Computational  space  C,  Parameter  space  P,  and  Physical  domain  D 

The  computational  space  C  is  defined  as  the  unit  square  in  a  two-dimensional  space  with  Carte¬ 
sian  coordinates  (£,77),  and  £  €  [0, 1],  7?  €  [0, 1]  (see  Figure  LI).  The  grids  are  uniformly  distributed 
on  the  boundaries  and  in  the  interior  area  of  the  computational  space.  The  mesh  sizes  are 
in  the  £  direction  and  'm  the  V  direction,  where  and  Nv  are  the  grid  numbers  in  the 
corresponding  direction,  ^he  parameter  space  V  is  defined  as  a  unit  space  in  a  two-dimensional 
space  with  Cartesian  coordinate  (s,  £),  and  s  €  [0, 1],  t  G  [0, 1].  The  boundary  values  of  s  and  t  are 
determined  by  the  grid  point  distribution  in  the  physical  domain. 


•  s  =  0  at  edge  Pi  and  s  ~  1  at  edge  P2 

•  s  is  the  normalized  arc-length  along  edges  P3  and  P4 


2 


20010221  137 


•  t  =  0  at  edge  £3  and  t  —  1  at  edge  E4 

•  t  is  the  normalized  arc-length  along  edges  E\  and  E? 

An  algebraic  transformation  s  :  C  — >  V  is  defined  to  map  the  computational  space  C  onto 
the  parameter  space  V.  The  grid  distribution  is  specified  by  this  algebraic  transformation,  which 
depends  on  the  prescribed  boundary  grid  point  distribution.  The  interior  grid  point  distribution 
inside  the  domain,  generated  by  the  algebraic  transformation,  is  a  good  reflection  of  the  prescribed 
boundary  grid  point  distribution.  Let  sg3(f)  =  s(£,  0)  and  S£4(£)  =  s(£,  1)  denote  the  normalized 
arc-length  along  edges  £3  and  £4j  tEl{y)  =  t{ 0,77)  and  tE2{y)  =  t(  1, 17)  denote  the  normalized 
arc- length  along  edges  E\  and  £2*  The  algebraic  transformation  s  :  C  — >  V  is  defined  as 


t  = 

Equation  (1.1)  is  called  the  algebraic  straight  line  transformation.  It  defines  a  differentiable  one- 
to-one  mapping  because  of  the  positiveness  of  the  Jacobian:  s^tv  —  >  0. 

The  elliptic  transformation  x  :  V  -*>  P,  which  is  independent  of  the  prescribed  boundary  grid 
point  distribution,  is  defined  to  map  the  parameter  space  V  onto  the  physical  domain  V.  The 
elliptic  transformation  is  equivalent  to  a  set  of  Laplace  equations 


Sxx  +  Syy  “  0  ^ 

i>XX  “f"  tyy  =  0 

The  elliptic  transformation  defined  by  the  above  equations  is  also  differentiable  and  one-to-one. 

Till  now  we  have  defined  two  transformations,  i.e.,  the  algebraic  transformation  s  :  C  —¥V,  and 
the  elliptic  transformation  x  :  V  -¥  V,  Because  both  the  algebraic  transformation  and  the  elliptic 
transformation  are  differentiable  and  one-to-one.  The  composition  the  two  transformation  is  also 
differentiable  and  one-to-one,  so  as  to  the  inverse  transformation. 

In  physical  domain,  the  curvilinear  coordinate  system  satisfies  a  system  of  Laplace  equations: 

Ar  =  0  (1.3) 

where  r  =  (x,y)T.  The  inherent  smoothness  of  the  Laplace  operator  makes  the  grids  smoothly 
distributed  in  the  physical  domain.  Being  transformed  to  the  computational  space,  this  Laplace 
system  becomes  a  set  of  Poisson  equations.  The  control  functions  is  determined  by  the  composed 
transformation  according  to  the  following  procedures. 

First,  Eq.(1.2)  is  transformed  into  the  computational  space  C: 


As  =  gnstf  +  2  g12szv  +  </22s^  +  Afsf  +  A  r)sv 
At  =  <7U%  +  2g12t£V  +  g22^  +  Aff^  +  Ar}tv 
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(1.4) 


where  <7U,<712,  J722  are  the  components  of  the  contravariant  metric  tensor,  which  can  be  calculated 
from  the  covariant  metric  tensor 

9n  =  j^922  =  {rv,r„)/J2 

9U  =  -^i2  =  -(r€,r^)/J2 

9 22  =  ^<7n  =  (r€,rf)/J2 

J  is  defined  as  ^det\gij\.  From  Eq.(1.2)  and  (1.4),  we  have 


and  the  matrix  T  is  defined  as 

T=(U”)  tu) 

Then  the  Laplace  system  Eq.(1.3)  is  transformed  to  the  computational  space  C: 

9nr((  +  2 gl2r(v  +  g72^  +  A£r(  +  A  r)rv  =  0.  (1.8) 

Substitute  Eq.(1.5)  into  Eq.(1.8),  A£  and  Ary  are  replaced  by  the  control  functions  on  the  right- 
hand-side  of  Eq.(1.5),  and  we  obtain  the  Poisson  equations  for  the  grid  generation  as  follows: 

,‘'rft+2S‘Vf, + +<r“*gVe + (s11^? +2<,‘!p<|> +9”pLV,  =  o 

(1.9) 

where  the  control  functions  ,  -^12^  ^22^ ?  -^11^ ^  -^12^  3  ^22^  3X6  determined  by  the  algebraic  trans¬ 
formation,  as  defined  previously  in  Eq.(1.6). 

The  elliptic  transformation  is  carried  by  solving  a  set  of  Poisson  equations.  The  control  functions 
are  specified  by  the  algebraic  transformation  only  and  it  is,  therefore,  not  needed  to  compute  the 
control  functions  at  the  boundary  and  to  interpolate  them  into  the  interior  of  the  domain,  as 
required  in  the  case  for  all  well-known  elliptic  grid  generation  systems  based  on  Poisson  systems. 

The  computed  grids  are  in  general  not  orthogonal  at  the  boundary.  The  algebraic  transformation 
can  be  redefined  to  obtain  a  grid  which  is  orthogonal  at  the  boundary.  First,  redefine  the  elliptic 
transformation  by  imposing  the  following  boundary  conditions  for  s  and  t: 
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•  s  —  0  at  edge  E\  and  a  —  1  at  edge  E^ 

•  =  0  along  edges  E3  and  #4,  where  n  is  the  outward  normal  direction 

•  t  =  0  at  edge  E3  and  t  =  1  at  edge  E$ 

i 

•  Ji  =  0  along  edges  E\  and  E2,  where  n  is  the  outward  normal  direction 

Second,  redefine  the  algebraic  transformation  s:C  ^  V  according  to  two  algebraic  equations, 


(1.10) 


(1.11) 


a  =  sEi(OH0(t)  +  sEi{OHi(t) 
t  =  tEl  {r))H0  (5)  +  tE 2  (a) 

where  Hq  and  H\  are  cubic  Hermite  interpolation  functions  defined  as 

H0(s)  =  (1  +  2s)(l  —  a)2 
ifi(a)  =  (3-2  a)a2 

Grid  orthogonahty  at  boundaries  is  obtained  in  three  steps. 

1.  Compute  an  initial  grid  based  on  the  Poission  grid  generation  system  with  control  functions 
specified  according  to  the  algebraic  straight  line  transformation  defined  by  Eq.1.1; 

2.  Solve  the  two  Laplace  equations  given  by  Eq.1.2  with  the  above  specified  boundary  conditions; 

3.  Recompute  the  g^lu  based  on  the  Poisson  system  but  with  control  functions  specified  according 
to  the  algebraic  transformation  defined  by  Eq.1.10. 

The  two-dimensional  grid  near  the  leading  edge  of  a  Joukowsky  airfoil  is  shown  in  Figure  1.2. 
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Figure  1.2:  Grid  near  the  leading  edge  of  a  Joukowsky  airfoil 
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1.2  Three-Dimensional  Grid  Generation 

The  basic  idea  of  three-dimensional  grid  generation  is  similar  to  that  of  the  two-dimensional  case. 
The  computational  space  is  a  unit  cubic  with  f  G  [0, 1],  17  G  [0, 1] ,  C  €  [0, 1].  The  parameter  space  is 
a  unit  cubic  with  s  G  [0,  l],t  G  [0, 1],»  G  [0, 1],  see  Figure  1.3.  « 

•  s  =  0  at  face  F\  and  s  =  1  at  face  £2 

•  3  is  the  normalized  arc-length  along  edges  Ei,  £2,  £3  and  £4 

•  t  =  0  at  face  £3  and  t  =  1  at  face  £4 

•  t  is  the  normalized  arc-length  along  edges  £5,  E$,  £7  and  £g 

•  u  —  0  at  face  £5  and  t  =  1  at  face  £6 

•  t  is  the  normalized  arc-length  along  edges  £9,  £10,  £11  and  £12 


Figure  1.3:  Computational  space  C,  Parameter  space  £,  and  Physical  domain  D 

Let  se^O  =  s(£,0,0),  SEj(f)  =  s(£,l,0),  se3(0  =  s(f,0,l),  and  sjs4(0  =  s(f,l,l)  denote 
the  normalized  arc-length  along  edges  E\,  £2,  £3,  and  £4;  t;;5 (r/)  =  t(0, r/,  0).  tE6{v )  =  t(l,r),0), 
tEhiv)  =  i(0, 77, 1),  and  ts8(j 7)  =  t(l,rj,  1)  denote  the  normalized  arc-length  along  edges  £5,  E6, 

£7  and  £8;  ttE„(0  =  **(0,0,0,  u£io(C)  =  *(l,0,O>  uEn(0  =  **(0,1,0)  and  ujs12(0  =  u(l»l>*?) 

denote  the  normalized  arc-length  along  edges  £9,  £10,  £11  and  £12.  The  algebraic  transformation 
from  computational  space  to  parameter  space  s  :  C  -t  V  is  defined  as 


s  =  sb,(0(1-^)(1-*1)  +  5Sj(^)*(1-*1)  +  sS3(6(1-*)**  +  «B4(^ 
t  =  tE6(0(l-«)(l-«)  +  ^6(,7)s(1-u)  +  ^(,?)(1-s)u  +  tJS8(,/)s*i  (L12) 

«  =  UE9(0(l-s)(l-*)  +  **Sio(0«(1-*)  +  u£:11(0(l-s)i  +  «E12(0^ 
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Equation  1.12  is  called  the  algebraic  straight  line  transformation. 

The  elliptic  transformation  x  :  V  -¥  P,  which  is  independent  of  the  prescribed  boundary  grid 
point  distribution,  is  defined  to  map  the  parameter  space  V  onto  the  physical  domain  V.  The 
elliptic  transformation  is  equivalent  to  a  set  of  Laplace  equations 


&XX  "b  Syy  +  SZZ 
txx  "b  tyy  ~b  ^ zz 

UXX  +  Uyy  +  Uzz 


0 

0 

0 


(1.13) 


The  elliptic  transformation  defined  by  the  above  equations  is  also  differentiable  and  one-to-one. 

In  physical  domain,  the  curvilinear  coordinate  system  satisfies  a  system  of  Laplace  equations: 

Ar  =  0  (1-14) 

where  r  =  (x,y,z)T. 

When  the  Laplace  system  Eq.(1.14)  is  transformed  to  the  computational  space  C,  we  obtain  a 
Poisson  system: 

9llr(f  +  J22^  4-  ff33rCf  +  2  g12riv  +  2  g13r^  +  2g23rf)C 
+  (9nPn  +  922P22  +  9™  Pm  +  l9UPn  +  2g13Ph  +  2</23P213)rf 

+  (9UP 11  +  9™Ph.  +  9^  Pm  +  2gUPh  +  2gl3P?3  + 

+  (9UP?i  +  9™P22  +  9™ Pm  +  2gUPn  +  2 g13P?3  +  2<723P233)rc  -  0  (1.15) 

where  ^l1,^22,^33,^12,^13,^23  are  contravariant  metric  tensor,  which  are  calculated  through  the 
co variant  metric  tensor 


9U  =  j2  (022033  -923) 


.22  _ 


.33 


=  ^(511033  -  513) 

(011022  -  9l2) 


912  =  J2(9l3923  -  912933) 

913  —  j2[gi2923  ~  913922) 

-  -^(ffl25l3  ~  02301l) 


r.23 


The  control  functions  are  defined  through  the  composite  transformation. 


P11  = 


P12  -  -T 


P22  =  -T 


Pi3  -  -T 


P33  =  -T-1 


P23  -  -T 


\  ) 


(1.16) 
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and  the  matrix  T  is  defined  as 


S£  3f] 
tri 


(1.17) 


The  computed  grids  are  in  general  not  orthogonal  at  the  boundary.  The  algebraic  transformation 
can  be  redefined  in  a  similar  way  as  in  the  previous  section  to  obtain  a  grid  which  is  orthogonal  at 
the  boundary.  The  three-dimensional  grid  around  a  delta  wing  is  shown  in  Figure  1.4. 


Figure  1.4:  Grid  around  a  delta  wing 
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Chapter  2 


Parallel  Computation 


2.1  Introduction 

Parallel  computing  provides  an  efficient  solution  for  large  scale  computations,  such  as  direct  nu¬ 
merical  simulation  (DNS)  and  large  eddy  simulation  (LES),  which  require  both  large  volume  of 
memory  and  high  speed.  In  a  parallel  computer,  several  processors  can  work  concurrently  on  vari¬ 
ous  parts  of  the  the  same  computation.  Some  parallel  computers,  such  as  SGI  Origin  2000,  have  a 
global  memory  that  can  be  accessed  by  all  the  processors.  Such  a  computer  is  said  to  have  a  share 
memory.  Most  programming  languages  and  environments  shield  the  programmer  from  working 
directly  with  processor.  Instead  they  supply  a  higher  lever  concept,  called  a  process,  that  models 
the  activation  and  a  single  program  on  a  processor.  The  shared-memory  programming  models 
are  characterized  by  the  provision  of  a  global  memory  that  can  be  directly  read  from  and  written 
to  by  every  process  involved  in  a  computation.  Compared  with  a  distributed-memory  system,  a 
shared-memory  computer  ran  provides  a  larger  band-width  for  communications  between  different 
processors. 

Two  families  of  programming  models  have  become  widely  used.  High  Performance  Fortran 
(HPF)  is  based  on  a  shared-memory  data-parallel  model.  The  other  major  parallel  programming  is 
a  distributed- memory  model  with  explicit  control  parallelism,  also  referred  to  as  a  message  passing 
programming  model.  Processes  in  a  message  passing  programming  model  are  only  able  to  read 
and  write  into  their  respective  local  memory.  They  synchronize  with  one  another  by  explicitly 
calling  library  procedures.  Data  is  copied  across  local  memories  by  having  the  appropriate  process 
send  and  receive  message  via  explicit  procedure  of  subroutine  calls.  The  Message  Passing  Interface 
(MPI)  standard  defines  a  set  of  functions  and  procedures  that  implements  the  message  passing 
model.  Although  the  message  passing  model  was  originally  designed  for  the  distributed-memory 
system,  it  can  also  be  used  for  a  shared-memory  system.  A  message  passing  model  built  on  a 
shared-memory  infrastructure  can  achieve  a  higher  speed  of  data  communication  because  of  the 
larger  band-width. 

In  a  shared-memory  multi-processor  system,  such  as  a  SGI  Origin  2000  computer,  the  speed  for 
a  processor  to  access  the  memory  on  its  node-board  (local  memory)  is  different  from  the  speed  for 
it  to  access  the  memory  of  other  processor  (remote  memory).  Therefore,  even  on  a  shared-memory 
system,  there  still  will  be  some  overhead  caused  by  data  exchange  between  different  processors. 
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2.2  Parallel  Scheme 


; 

In  our  computation,  the  parallel  computing  is  based  on  MPI.  The  computer  domain  is  divided  into 
Np  sub-domains  along  the  £  direction,  where  Np  is  the  number  of  processes.  The  data  of  each  sub- 
domain  is  stored  within  the  local  memory  of  each  po}ce$sor.  Each  processor  communicates  with 
each  other  using  a  group  functions  provided  by  MPT  library.  Data  exchange  between  processors 
occurs  when  the  data  stored  in  the  local  memory  of  other  processor  is  required.  For  instance, 
in  order  to  calculate  the  first  and  second  derivative  in  the  £  direction  using  the  compact  scheme, 
massive  data  exchange  is  achieved  by  calling  the  MPI  function  MPI_ALLTOALL.  In  some  other 
case,  when  data  exchange  is  required  between  the  interface  of  two  adjacent  sub-domain,  the  MPI 
function  MPIJ3END  and  MPLRECV  are  used. 


2.3  Parallel  Speedup 


A  test  code  has  been  used  to  compare  the  performance  of  parallel  computing  on  a  SGI  Origin  2000 
and  a  Cray  T3E.  The  SGI  Origin  2000  is  a  single  large  shared-memory  system  which  is  also  capable 
to  run  the  MPI  code  designed  for  a  multiple-instruction  multiple-data  (MIMD)  system.  The  SGI 
Origin  2000  we  are  testing  has  sixteen  R10000  MIPS  processors  at  250MHz  and  4GB  of  memory. 

A  MPI  code  is  designed  to  test  the  performance  of  computers.  In  this  testing  code,  first 
derivative  in  £,  7]  and  £  directions  is  computed  using  the  4th  order  compact  scheme,  which  is 
popular  high-order  scheme  in  DNS/LES. 

The  grid  numbers  in  the  three  direction  are  480,  160,  and  80  respectively,  which  represent  a 
moderate  scale  computation.  Several  parameters  are  used  to  measure  the  performance  of  a  parallel 
computing.  The  wall-clock  time  is  the  elapsed  time  of  a  MPI  program  which  can  be  obtained  by 
calling  MPI  subroutine  MPI.WTIME.  Speedup  is  Penned  as  the  ratio  of  the  runtime  of  a  serial 
solution  to  a  problem  to  the  parallel  runtime. 


S{n,p) 


Tsjn) 

T„{n,p) 


(2.1) 


where  Tj(n)  denotes  the  runtime  of  the  serial  program  running  with  one  process.  Tn(n,p)  denotes 
the  runtime  of  the  parallel  code  running  with  p  processes. 


•  S(n,p)  >  p:  the  parallel  program  is  exhibiting  super- linear  speedup; 

•  S(n,p)  =  p:  the  parallel  program  is  exhibiting  linear  speedup; 

•  1  <  S(n,p)  <  p:  the  parallel  program  is  exhibiting  speedup; 

•  5(n,p)  <  1:  the  parallel  program  is  exhibiting  slowdown; 


An  alternative  to  speedup  is  efficiency,  which  is  a  measurement  of  process  utilization  in  a  parallel 
program,  relative  to  the  serial  program.  It  is  defined  as 


E{n,p)  = 


S(»,P) 

P 


Tf{n) 

pTv{n,p) 


(2.2) 
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•  E{n,p)  >  1:  the  parallel  program  is  exhibiting  super-linear  speedup; 

•  E(n,p)  =  1:  the  parallel  program  is  exhibiting  linear  speedup; 

•  1/p  <  E(n,p)  <  1:  the  parallel  program  is  exhibiting  speedup; 

•  E(n,p)  <  1/p:  the  parallel  program  is  exhibiting  slowdown; 

Figure  2.1  presents  the  wall  clock  time  for  2,  4,  8  16,  and  32  processors  on  Cray  T3E.  Figure 
2.2  shows  a  similar  result  obtained  on  the  SGI  Origin  2000  computer.  On  both  machines,  the  wall 
clock  decreases  as  the  number  of  processors  increases.  The  comparison  is  displayed  in  Figure  2.3. 
The  wall  clock  time  to  run  the  same  case  on  T3E  is  about  2/3  of  the  time  on  the  Origin  2000  for 
2,  4,  8,  and  16  processors. 


Figure  2.1:  Wall-clock  time  for  different  num-  Figure  2.2:  Wall-clock  time  for  different  num¬ 
ber  of  processors  on  Cray  T3E  ber  of  processors  on  SGI  Origin  2000 


Figure  2.3:  Wall-clock  time  for  different  number  of  processors  on  Cray  T3E  and  SGI  Origin  2000 

The  performances  indicated  by  speedup  of  MPI  code  running  on  T3E  and  Origin  2000  are 
displayed  in  Figure  2.4  and  2.5.  The  computing  of  derivative  in  r]  and  £  direction  display  super- 
linear  scalability.  Because  there  is  no  data  exchange  in  rj  and  £  direction,  the  parallel  computing 
is  very  efficient  to  calculate  derivative  in  these  two  directions.  On  Origin  2000,  the  super-linear 
speedup  can  be  explained  as  following.  As  the  number  of  processor  increase,  the  ratio  between  the 
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cache  of  each  processor  and  the  data  handled  by  the  processor  also  increase,  in  other  words,  more 
data  can  be  put  into  the  cache,  which  has  a  faster  accessing  speed. 

But  on  T3E,  only  sub-linear  result  is  obtained.  When  we  submitted  the  testing  MPI  job  on  SGI 
Origin  2000,  no  other  job  was  running.  But  our  MPI  testing  code  was  submitted  to  T3E  running 
with  many  jobs  submitted  by  other  users,  and  only  wall-clock  time  has  been  recorded.  The  speedup 
parameter  does  not  display  the  linear  scalability  for  calculating  derivatives  in  £  direction,  because 
massive  data  exchange  is  required  by  calling  the  MPI  function  MPI-A.LLTOALL.  The  data  exchange 
between  different  processors  has  reduced  the  totcil  performance  of  parallel  computing.  This  extra 
cost  for  parallel  computing  is  inevitable.  The  speedup  curves  as  a  function  of  number  of  processors 
for  T3E  and  Origin  are  displayed  in  Figure  2.6.  The  total  speedup  of  T3E  is  better  than  that  of 
Origin  2000. 


Figure  2.4:  Speedup  for  different  number  of  Figure  2.5:  Speedup  for  different  number  of 
processors  on  Cray  T3E  processors  on  SGI  Origin  2000 


Figure  2.6:  Speedup  for  different  number  of  processors  on  Cray  T3E  and  SGI  Origin  2000 

The  efficiency  parameter  are  displayed  in  Figure  2.7  and  2.8.  On  Cray  T3E,  when  the  number 
of  processors  is  less  than  4,  the  efficiency  for  calculating  derivative  in  £  direction  and  the  total 
performance  is  greater  than  1.  That  means  even  the  T3E  is  loaded  with  jobs  submitted  from 
different  users,  four  processors  are  always  available  for  our  job.  Therefore,  super  linear  performance 
is  obtained.  On  Origin  2000,  no  matter  how  many  processors  are  used,  the  efficiency  for  calculating 


13 


E(".P) 


derivative  in  r)  and  (  direction  is  always  greater  than  1,  which  means  super  linear  performance 
attained. 


Number  of  processors  Number  of  processors 


Figure  2.7:  Efficiency  for  different  number  of 
processors  on  Cray  T3E 


Figure  2.8:  Efficiency  for  different  number  of 
processors  on  SGI  Origin  2000 


Chapter  3 

Numerical  Simulation  of  Separated  Flow  around  a 
NACA  0012  airfoil  at  12°  Angle  of  Attack 


3.1  Introduction 

It  is  of  particular  interest  to  study  the  flow  separations  around  airfoils  at  large  angle  of  attack.  Flow 
separations  have  at  least  two  effects.  The  first  one  is  the  sudden  loss  of  lift,  and  the  second  one  is 
the  generation  of  aerodynamic  noise.  These  two  aspects  are  crucial  to  the  design  of  aircraft,  and 
these  problems  can  not  be  solved  in  the  absence  of  the  detailed  information  about  flow  separation, 
which  is  a  complex  time-dependent  physical  process.  Understanding  of  the  process  is  still  an  open 
question  for  research.  However,  this  kind  of  problem  has  importance  in  instantaneous  quantities 
and  must  studied  by  DNS  or  LES  at  least,  but  not  by  RANS  (Reynolds-averaged  Navier  Stokes). 

The  spatial  and  temporal  complexity  of  this  problem  makes  it  inaccessible  by  conventional 
experimental  and  numerical  techniques.  From  the  experimental  point  of  view,  as  it  has  been 
pointed  out  by  Shih  et  al.  (1992),  that  the  level  of  understanding  needs  to  advance  from  qualitative 
conjectures  based  on  now  visualization  and/or  measurement  of  global  quantities,  such  as  lift  and 
drag,  to  the  quantitative  measurement  of  the  instantaneous  flow  field.  The  prominent  feature 
about  the  flow  past  an  airfoil  at  large  angles  of  attack  is  the  emergence  of  large-scale  vortices 
when  the  flow  separates  from  the  airfoil  surface.  The  spatial  and  temporal  evolution  of  these 
vortical  structures  dominates  the  unsteady  flow  behavior  over  the  airfoil.  Detailed  information, 
such  as  spatial  vorticity  distribution  at  each  instance,  is  required  to  study  in  this  problem.  This 
requirement  excludes  the  use  of  traditional  single-point  velocity  measurement  techniques,  such  as 
hot-wire  anemometry  or  laser  Doppler  anemometry.  A  new  experimental  technique,  particle  image 
displacement  velocimetry  (PIDV),  has  been  used  by  Shih  et  al.  (1992)  to  study  the  unsteady  flow 
past  a  NACA  0012  airfoil  in  pitching-up  motion  in  a  water  towing  tank.  The  Reynolds  number 
based  on  the  airfoil  chord  and  the  freestream  velocity  is  5000.  Instantaneous  velocity  fields  at 
different  time  have  been  acquired  over  the  entire  2-D  flow  field.  Using  this  flow  field  data,  the 
out-of-plane  (spanwise)  component  of  vorticity  is  computed  to  show  the  vortical  structures.  They 
found  out  that  boundary-layer  separation  near  the  airfoil  leading-edge  leads  to  the  formation  of  a 
vortical  structure.  Near  the  leading-edge,  vorticity  accumulations  occurs  because  of  the  slowdown 
of  the  downstream  convection  as  a  result  of  the  adverse  pressure  gradient.  The  accumulation 
process  is  eventually  interrupted  by  a  sudden  emergence  of  unsteady  flow  separation,  with  the 
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immediate  release  of  the  accumulated  vorticity  into  the  outer  flow.  This  process  may  be  explained 
by  the  separation  model  proposed  by  Van  Dommelen  fc  Shen  (1980)  and  Van  Dommelen  &;  Cowley 
(1990).  This  model  describes  the  process  of  initial  breakup  of  unsteady  boundary  layer  through 
a  significant  thickening  of  the  boundary  layer.  In  this  case,  the  thickening  of  boundary  layer  is 
driven  by  the  reversing  fluid  layer  which  is  formed  as  a  result  of  the  adverse  pressure  gradient. 
In  the  reversed  flow  region,  a  streamwise  accumulation  of  boundary  layer  vorticity  occurs.  The 
eruption  of  the  boundary  layer  vorticity  lead  to  formation  and  shedding  of  large  scale  vortical 
structures  which  are  named  as  the  primary  vortices  from  the  leading-edge.  The  primary  vortices 
move  downstream  along  the  airfoil  surface  and  interact  with  the  local  boundary  layer  below  to  cause 
the  upstream  accumulation  and^  eruption  of  reversal  boundary  layer  vorticity.  Near  the  trailing- 
edge,  the  primary  vortex  further  triggers  the  formation  and  shedding  of  a  count-rotating  vortex. 
The  unsteady  separated  flow  near  the  leading-edge  and  trailing  edge  of  a  NACA  0012  airfoil  was 
measured  in  more  details  by  Shih  et  al.  (1995)  using  particle  image  velocimetry  (PIV)  technique. 
In  his  experiment,  the  Reynolds  number  was  5000  and  25000,  based  on  the  chord  length  of  the 
airfoil  and  the  freestream  velocity.  The  role  of  the  trailing  edge  flow  was  also  examined.  At  the  late 
stage,  as  the  primary  vortex  approach  the  trailing  edge,  a  counter-rotating  vortex  is  shed  from  the 
trailing  edge.  The  experimental  results  obtained  by  Shih  et  al  (1992,  1995)  will  be  compared  with 
our  DNS  results  in  Section  3.  It  is  believed  by  the  authors  that  the  generic  characteristics  of  an 
unsteady  separated  flow  are  fairly  universal  and  independent  of  the  Reynolds  number  and  external 
flow  conditions.  The  Reynolds  number  has  effect  on  the  time  and  length  scales  of  the  separation 
structure.  For  example,  the  separation  structure  of  a  low  Reynolds  number  flow  evolves  faster  than 
a  high  Reynolds  number  case.  The  stronger  viscous  diffusion  of  a  low  Reynolds  number  flow  will 
thicken  the  boundary  layer  and  thereby  attenuate  the  explosive  nature  of  the  unsteady  separation 
(Shih  et  al.  1995). 

For  numerical  approaches,  Tenaud  &  Phuoc  (1997)  used  large  eddy  simulation  (LES)  to  study 
separated  flow  around  a  NACA  0012  airfoil  at  20°  angle  of  attack.  Three  different  flow  regions 
following  different  structure  behaviors  were  obse^ed  in  their  results.  The  three  regions  are  the 
leading-edge,  the  middle  part  of  the  upper  surface,  and  the  trailing-edge.  The  leading-edge  region 
is  dominated  by  vortex  shedding  due  to  separation  of  the  boundary  layer.  The  second  region 
corresponds  to  the  middle  part  of  the  upper  surface  of  the  airfoil  where  the  eddy  structures  grow 
and  are  convected  downstream.  The  last  region  is  situated  close  to  the  trailing-edge  where  alternate 
vortices  are  created  due  to  the  wake  instability.  The  authors  did  not  provide  a  detailed  description 
of  the  procedure  of  vortex  shedding  from  the  leading-edge  and  the  subsequent  evolution  of  vortical 
structures. 

This  work  focuses  on  direct  numerical  simulation  of  flow  separation  around  a  NACA  0012 
airfoil  at  a  12°  angle  of  attack.  In  this  case,  the  following  understanding  has  been  obtained  from 
the  previous  theoretical  and  experimental  investigation.  First,  the  large  vortices  are  intermittently 
formed  and  shed  from  the  leading-edge.  Second,  the  separation  leads  to  an  alternation  of  the 
pressure  distribution  and  therefore  changes  the  lift  and  moment  acting  on  an  airfoil.  Thus,  the 
aerodynamic  forces  become  unsteady  and  there  is  a  dramatic  decrease  in  lift  accompanied  by  an 
increase  in  drag  and  large  changes  in  the  moment  exerted  on  the  airfoil.  Third,  the  interactions 
between  vortices  and  between  vortices  and  airfoil  surface  will  generate  noise.  However,  a  complete 
understanding  of  these  complex  unsteady  effects  has  not  yet  been  achieved,  and  there  is  a  great 
need  for  systematic  fundamental  studies.  Understanding  of  flow  separation  will  provide  assistance 
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to  improving  the  separation  control  and  noise  control.  The  main  objectives  of  the  present  work  is 
to  reveal  the  detailed  flow  separation  structures  using  high-resolution  simulation. 


3.2  Governing  Equations  And  Numerical  Methods 


The  two-dimensional  Navier-Stokes  equations  in  generalized  curvilinear  coordinates  (£,  rj)  are  writ¬ 
ten  in  conservative  forms: 

i  nn  H(  z 7  _  p.t  \  HI  h'  —  \ 
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0  \ 
Txx^x  “b  ryx£y 
TXy£x  +  Tyy£y 


Fv  =  T 


TxxVx  H“  TyxVy 
TxyVx  TyyVy 
\  qxVx  +  QyVv  } 


where  J  is  Jacobian  of  the  coordinate  transformation,  and  £x,€y>Vx,i Ty  are  coordinate  transfor¬ 
mation  metrics,  p  is  density,  p  is  pressure,  tx  and  v  are  components  of  velocity.  U  =  u£x  +  v(y, 
V  —  UTjx  +  VTjy.  e  is  the  total  energy.  The  components  of  viscous  stress  and  heat  flux  are  denoted 
by  Xu,  Tjy,  Tjiy  and  (jx ,  Qy ,  respectively. 

In  Eq.  (3.1),  the  second  order  Euler  Backward  scheme  is  used  for  time  discretization,  and  the 
fully  implicit  form  of  the  discretized  equation  is: 


3Qn+1  -  4 Qn  +  QnJ^_  d{En+\  -  E?+1) 


+ 


2JAt 
8{Fn+ 1  -  F£+1) 

dv 


dS 


=  0 


(3.2) 


Qn+1  is  estimated  iteratively  as: 


Qn+1  =  Qp  +  SQP 

SQP  =  QP+ 1  -  Qp 

At  step  p  =  0,  QP  =  Qn]  as  5QF  is  driven  to  zero,  Qp  approaches  Qn+1.  Flux  vectors  axe  linearized 
as  following: 

E"+1  *  EP  +  APSQP 
Fn+1  «F  +  Bp5Qp 
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(3.3) 


So  that  Eq.  (3.2)  can  be  written  as: 

[§/  +  A  tJ{D^A  +  DVB)]SQP  =  R 

where  R  is  the  residual: 

R  =  -  2 Qn  +  ig"-1)  -  A tJ[{D({E  -  Ev)  +  DV(F  -  F„)f 

The  superscript  p  stands  for  iteration  step.  D £,  D ^  represent  partial  differential  operators,  A,  B 
are  the  Jacobian  matrices  of  flux  vectors: 


The  right  hand  side  of  Eq.  (3.3)  is  discretized  using  the  fourth-order  compact  scheme  (Lele, 
1992)  for  spatial  derivatives,  and  the  left  hand  side  of  the  equation  is  discretized  following  LU-SGS 
method  (Yoon  &  Kwak,  1992).  In  this  method,  the  Jacobian  matrices  of  flux  vectors  are  split  as: 

A  =  +  B  =  B++B~ 


where, 

A*  =  \[A±taI\ 

B±  =  \[B±tbI\ 

where, 

rA  =  «maa:[|A(A)|]  +  v 
tb  =  nmax[\X(B)\)  +  u 

where  A(A),  \{B)  are  eigenvalues  of  A,B  respectively,  n  is  a  constant  greater  than  1.  V  represents 
the  effects  of  viscous  terms.  The  following  expression  is  used. 

"  =  ma^{1-l)M?RePT’  (3'4) 


The  first-order  upwind  finite  difference  scheme  is  used  for  the  split  flux  terms  on  the  left  hand 
side  of  Eq.  (3.3).  This  does  not  effect  the  accuracy  of  the  final  solution.  As  the  left  hand  side  is 
driven  to  zero,  the  discretization  error  of  the  left  hand  side  will  also  be  driven  to  zero'  The  finite 
difference  representation  of  Eq.  (3.3)  can  be  written  as: 

+  &tJ{rA  +  rB)I]6(^ij  = 

RFi.-teJ[A-6Q?i+lj-A+6Ct>i_ld 

+BS%+l-B+S%_ i  ] 

In  LU-SGS  scheme,  Eq.  (3.5)  is  solved  by  three  steps.  First  we  initialize  6Q°  using 

8Ql  =  [\l  +  MJ{rA  +  TB)I]-1^  '  '  (3.5) 
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In  the  second  step,  the  following  formula  is  used: 


so*  =  SQh 

+  [|j  +  A  tJ(rA  +  rB)/]-1 

x  [AUiA+SQUj  +  B+SQtj^)]  (3.6) 

For  the  last  step,  SQP  is  obtained  by 

6%  = 

-  [|  I  +  AtJ(rA  +  rB)]-1 

x  [A  tJ{A-8^+B-6^+x))  (3.7) 

The  sweeping  is  performed  along  the  planes  of  i  +  j  =  const ,  i.e  in  the  second  step,  sweeping  is 
from  the  low-left  corner  of  the  grid  to  the  high-right  comer,  and  then  vice  versa  in  the  third  step. 

In  order  to  depress  numerical  oscillation  caused  by  central  difference  scheme,  spatial  filtering 
is  used  instead  of  artificial  dissipation.  The  implicit  sixth-order  compact  scheme  for  space  filtering 
(Lele,  1992)  is  applied  for  primitive  variables  u,  u,p,p  afte»  each  time  step. 

For  subsonic  flow,  u,  v,  T  are  prescribed  at  the  upstream  boundary,  p  is  obtained  by  solving  the 
modified  N-S  equation  based  on  characteristic  analysis.  On  the  far  field  and  out  flow  boundary, 
non-reflecting  boundary  conditions  are  applied.  Adiabatic,  non-slipping  condition  is  used  for  the 
wall  boundary.  All  equations  of  boundary  conditions  are  solved  implicitly  with  internal  points. 
Specific  details  of  boundary  treatment  can  be  found  in  Jiang  et  al .  (1999). 


3.3  Computational  Details 

Direct  numerical  simulation  has  been  implemented  to  investigate  the  compressible  flow  separation 
around  a  NACA  0012  airfoil  at  a  12°  angle  of  attack.  The  chord  length  L  is  taken  as  the  charac¬ 
teristic  length,  and  the  freestream  velocity  is  the  characteristic  velocity.  The  Reynolds  number 
based  on  the  chord  length  and  the  freestream  velocity  is  Rec  =  5  x  105.  The  freestream  Mach 
number  is  0.4. 

An  elliptic  grid  generation  method  first  proposed  by  Spekreijse  (1995)  is  used  to  generate  the 
2-D  C-type  grids.  The  elliptic  grid  generation  method  is  based  on  a  composite  mapping,  which 
consists  of  a  nonlinear  transfinite  algebraic  transformation  and  an  elliptic  transformation.  The 
algebraic  transformation  maps  the  computational  space  onto  a  parameter  space,  and  the  elliptic 
transformation  maps  the  parameter  space  on  to  the  physical  domain.  The  elliptic  transformation 
is  carried  by  solving  a  set  of  Poisson  equations.  The  control  functions  are  specified  by  the  algebraic 
transformation  only  and  it  is,  therefore,  not  needed  to  compute  the  control  functions  at  the  bound¬ 
ary  and  to  interpolate  them  into  the  interior  of  the  domain,  as  required  by  the  well-known  elliptic 
grid  generation  systems  based  on  Poisson  systems  (Thompson  et  al.7  1985).  The  orthogonality  of 
grids  on  the  body  surface  or  near  the  boundary  is  attainable  by  re-configuration  of  the  algebraic 
transformation. 
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In  Figure  3.1(a),  the  body-fitted  C-grid  around  a  NACA  0012  airfoil  is  displayed.  The  mesh 
near  the  airfoil  surface  is  shown  in  Figure  3.1(b)  where  grids  axe  orthogonal  at  the  boundaries.  The 
mesh  is  tightened  near  the  wall  in  the  wall-normal  direction  (77)  as  well  as  in  the  vicinity  of  the 
leading-  and  trailing-edge  of  the  airfoil  in  the  direction  parallel  to  the  wall  (£).  The  mesh  size  is 
841  points  in  the  £  direction  and  141  points  in  the  77  direction. 


(a)  overview  of  the  grids 


(b)  grid  near  airfoil  surface 


Figure  3.1:  C-grid  around  a  NACA  0012  airfoil 


3.4  Results  And  Discussions 

Flow  separations  around  an  airfoil  with  a  high  angle  of  attack  (a  =  12°)  has  been  studied  using 
high-resolution  numerical  simulation.  In  this  case,  the  fluid  flow  around  the  airfoil  becomes  very 
unstable  and  different  eddy  structures  are  formed  in  the  vicinity  of  the  airfoil.  These  eddy  structures 
severely  affect  the  airfoil’s  aerodynamic  properties.  Vortex  interactions  and  the  interactions  between 
the  airfoil  and  vortex  may  lead  to  noise  generation.  In  the  present  work,  effort  is  concentrated  on 
providing  a  detailed  picture  of  flow  separation  process. 

The  computation  starts  from  freestream  conditions  and  without  any  perturbation.  During  the 
numerical  simulation,  the  flow  separation  process  is  found  to  be  unsteady  and  the  flow  field  is 
recorded  every  1000  time  steps,  where  the  time  step  is  approximately  At  =  1.767  x  10 ~4L/Uoo, 
An  animation  of  the  flow  field  has  been  made  based  on  the  recorded  data  to  show  the  contours  of 
the  instantaneous  spanwise  vorticity.  Eighteen  snap-shorts  with  equally  spaced  time  intervals  sure 
extracted  from  the  animation  and  displayed  in  Figure  3.2.  The  time  interval  between  two  adjacent 
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frames  is  about  0.88AL/U.  For  each  frame  the  region  in  the  vicinity  of  the  trailing-edge  are  enlarged 
and  displayed  in  Figure  3.3.  In  both  Figure  3.2  and  Figure  3.3,  the  contours  of  positive  vorticity 
are  plotted  in  solid  lines  and  negative  vorticity  are  in  dotted  lines. 

In  frame  (a)  of  Figure  3.2,  instability  is  observed  only  in  the  wake  near  the  trailing-edge  of  the 
airfoil.  The  vicinity  of  trailing-edge  is  enlarged  and  shown  in  frame  (a)  of  Figure  3.3.  Due  to  the 
12°  angle  of  attack,  a  strong  shear  layer  forms  at  Ihe  trailing  edge  of  the  airfoil.  The  presence  of 
alternate  vortices  indicates  the  Kelvin-Helmholtz  type  instability,  which  is  well-known  in  the  wake 
where  free  shear  layer  appears. 

In  frame  (b)  of  Figure  3.2,  a  separation  bubble  can  be  observed  on  the  upper  surface  of  the  airfoil 
near  the  leading-edge.  The  leading-edge  separation  bubble  is  formed  when  the  laminar  boundary 
layer  separates  from  the  surface  as  a  result  of  the  strong  adverse  pressure  gradient  downstream  of  the 
point  of  minimum  pressure.  According  to  the  separation  model  proposed  by  Van  Dommelen  &  Shen 
(1980)  and  Van  Dommelen  &  Cowley  (1990),  the  boundary  layer  is  thickened  around  the  position  of 
vorticity  reversal  due  to  streamwise  accumulation  of  fluid  driven  by  the  adverse  pressure  gradient. 
The  accumulation  of  vorticity  eventually  leads  to  the  eruption  of  the  boundary  layer  vorticity  and 
shedding  of  large  scale  vortical  structures  which  are  named  as  the  primary  vortices.  The  shedding 
of  the  primary  vortices  can  be  observed  in  every  picture  after  frame  (c)  of  of  Figure  3.2.  These 
vortical  structures  are  rotating  in  the  clockwise  direction  and  with  positive  spanwise  vorticity.  A 
typical  case  of  leading-edge  vortex  shedding  is  presented  by  frame  (d)  of  Figure  3.2  where  five 
vortices  are  sequently  shedding  from  the  leading-edge  and  moving  downstream  due  to  convection. 
Similar  phenomena  are  captures  by  the  snapshot  in  frames  (i),  (1),  (m),  (q)  and  (r)  of  Figure  3.2. 
Below  the  primary  vortex  near  the  wall,  a  layer  of  reversed  vorticity  is  generated  due  to  the  induced 
motion  of  the  primary  vortex.  As  a  result,  the  primary  vortices  are  propelled  away  from  the  wall 
by  the  reversed  vorticity  near  the  surface. 

In  frame  (e)  of  Figure  3.2,  two  clockwise  rotating  primary  vortices  roll  around  each  other  at  the 
middle  of  the  airfoil’s  upper  surface  and  eventually  merge  into  one  large  vortex  with  the  original 
rotating  direction.  The  vortex  merging  is  assorted  with  the  adverse  pressure  gradient  along 
the  upper  surface  of  the  airfoil.  After  departure  from  the  leading-edge,  the  streamwise  propagation 
velocity  of  the  primary  vortices  is  decelerated  as  a  result  of  the  adverse  pressure  gradient.  Therefore, 
a  downstream  primary  vortex,  which  leaves  the  leading-edge  earlier  and  slows  down  due  to  the 
adverse  pressure  gradient,  will  be  caught  up  by  an  upstream  vortex  which  has  a  larger  convective 
speed.  Similar  process  are  seen  in  frames  (h),  (o),  and  (p)  of  Figure  3.2. 

The  vortex  originated  from  the  the  merging  of  two  primary  vortices  seems  to  re-attach  on  the 
upper  surface  of  the  airfoil.  It  interacts  with  the  local  boundary-layer  and  induces  an  upstream 
accumulation  and  eruption  of  a  reversed  boundary-layer  vorticity  from  the  wall,  see  frame  (f)  of 
Figure  3.2.  In  this  frame  another  event  of  vortex  merging  is  observed  near  the  one-third  chord 
position.  The  accumulated  reversed  vorticity  leaves  the  wall  as  a  result  of  the  induced  motion  by 
the  downstream  primary  vortex  and  evolves  into  a  vortex  with  negative  spanwise  vorticity,  shown 
in  dotted  isovorticity  contours  in  frame  (f)  of  Figure  3.2.  This  induced  vortex  interacts  with  the 
primary  vortex  to  form  a  counter-rotating  vortex  pair.  A  typical  vortex  pairing  phenomenon  is 
presented  in  frames  (g)  of  Figure  3.2  where  two  vortex  pairs  are  observed.  Similar  phenomena 
of  vortex  pairing  have  been  described  in  experiment  of  Shih  et  al.  (1992).  The  vortex  pairs  are 
propelled  from  the  wall  by  a  self-induced  motion,  whereas  they  encounter  a  faster  convective  motion 
and  move  downstream  more  quickly.  At  the  same  time,  the  vortex  pairs  experience  substantial 
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deformation  due  to  the  difference  in  streamwise  velocity  of  the  main  flow  close  to  wall  and  that  far 
from  the  wall.  Because  the  far  end  of  the  vortex  pair  moves  faster  than  the  near- wall  end,  the  vortex 
pair  is  stretched  in  the  streamwise  direction,  as  shown  in  frame  (h)  of  Figure  3.2.  The  vortex  pair 
is  dominated  by  the  clockwise  rotating  primary  vortex  which  is  stronger  than  the  reversed  vortex. 
And  streamwise  stretching  occurs  to  the  clockwise  rotating  primary  vortex.  In  frame  (h),  the 
downstream  vortex  pair  shown  in  frame  (g)  merges  with  an  another  vortex  pair  which  comes  from 
the  wall.  This  vortex  merging  can  be  seen  more  clearly  in  frame  (h)  of  Figure  3.3.  In  frame  (i) 
of  Figure  3.2,  ail  these  vortex  pairs  shown  in  frame  (g)  and  (h)  seem  to  merge  with  each  other  to 
form  a  more  complex  vortex  system  near  the  trailing  edge  of  the  airfoil. 

When  the  dominating  clockwise-rotating  part  of  the  vortex  system  sweeps  over  the  trailing 
edge  of  the  airfoil,  a  counterclockwise-rotating  vortex  which  can  be  seen  in  frame  (h)  of  Figure  3.3 
appears  at  the  trailing-edge.  This  vortex  is  named  as  the  trailing-edge  vortex,  which  is  generated 
due  to  the  induced  motion  of  the  clockwise-rotating  primary  vortex.  As  clockwise-rotating  primary 
vortex  reaches  the  trailing-edge,  it  introduces  a  local  low  pressure  region  which  drives  the  lower 
surface  of  the  separating  shear  layer  to  curve  upward  and  move  across  the  wake  into  the  upper 
stream.  Then,  induced  by  the  clockwise-rotating  primary  vortex,  the  shear  layer  quickly  rolls  into 
an  intense  counterclockwise-rotating  trailing-edge  vortex.  This  process  is  observed  clearly  in  frame 
(h)  and  (i)  of  Figure  3.3.  Similar  phenomenon  was  also  reported  in  experiments  (Shih  et  al.  1992, 
1995). 
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In  frame  (j)  of  Figure  3.3,  the  trailing-edge  vortex  is  pushed  away  from  the  wall  due  to  a  self- 
induced  motion  as  the  vortex  grows.  The  trailing-edge  shear  layer  seems  to  stick  to  the  downstream 
edge  of  the  large-scale  trailing-edge  vortex  and  is  stretched  as  the  trailing-edge  vortex  is  propelled 
away  from  the  wall.  A  chain  of  small-scale  counterclockwise-rotating  vortical  structures  appears 
along  the  shear  layer  due  to  the  Kelvin-Helmholtz  type  instability.  Finally  the  trailing-edge  vortex 
separates  from  the  shear  layer  to  form  a  vortex  pair  with  a  clockwise-rotating  primary  vortex  and 
diffuses  away  from  the  trailing-edge,  as  shown  in  frame  (j)  and  (k)  of  Figure  3.2.  As  the  trailing- 
edge  vortex  leaves  the  trailing  edge,  pressure  recovers  near  upper  surface.  Therefore,  the  shear 
layer  moves  back  and  settles  along  the  streamwise  direction,  see  frames  (k)  and  (1)  of  Figure  3.3. 
The  Kelvin-Helmholtz  type  instability  structures  are  still  visible  along  the  shear  layer  until  another 
large-scale  primary  vortex  approaches  the  trailing-edge  and  leads  to  the  generation  of  another 
trailing-edge  vortex. 

In  the  results  presented  in  this  paper,  the  leading-edge  vortex  continues  to  shed  and  convect 
downstream.  Along  the  upper  surface  of  the  airfoil,  vortex  merging  occurs  due  to  the  adverse 
pressure  gradient.  The  leading-edge  primary  vortex  interacts  with  local  boundary  layer  and  triggers 
reversed  vorticity  accumulation  and  eruption,  which  in  turn  evolves  into  the  counterclockwise- 
rotating  vortical  structures  to  form  vortex  pairs  with  the  primary  vortex.  When  these  large-scale 
vortical  structures  arrive  the  trading-edge,  a  reversed  vortical  structure  is  induced  from  the  trailing- 
edge.  This  counterclockwise-rotating  vortex  moves  ' away  from  the  wall  by  a  self-induced  motion 
and  separate  with  the  shear  layer  to  interact  with  the  primary  vortex  and  form  a  vortex  pair, 
which  diffuses  away  from  the  trailing-edge  and  moves  downstream.  Therefore,  the  generation  of 
the  trailing-edge  vortex  is  a  quasi-periodic  process  which  can  be  seen  from  frame  (a)  through  (r) 
of  Figure  3.2  and  Figure  3.3. 

To  obtain  the  spectrum  information  about  flow  separation,  twelve  points  are  selected  near  the 
airfoil  surface  and  in  the  wake  to  record  the  time  series  of  instantaneous  velocity  and  pressure.  The 
location  of  these  points  are  denoted  by  PI  -  P12  and  shown  in  Figure  3.4.  Points  PI  -  P10  locate 
above  and  very  close  to  the  suction  surface  of  the  airfoil.  Points  Pll,  P12  are  in  the  wake. 


Figure  3.4:  Location  points  PI  -  P12  where  the  time  series  are  recorded 

The  instantaneous  fluctuations  of  the  streamwise  velocity  and  pressure  at  location  point  PI  and 
their  power  spectrum  are  displayed  in  Figure  3.5  -  Figure  3.8.  Here,  the  mean  flow  are  defined  based 
on  the  temporal  average.  Figure  3.6  shows  the  power  spectrum  of  streamwise  velocity  fluctuations, 
where  the  first  peak  associated  with  the  slow  variation  of  mean  flow  is  ignored.  The  frequency  of 
the  second  and  the  third  peak  is  /  =  1.28,  /  =  1.81.  In  this  case,  the  vortices  shed  from  the  leading 
edge  intermittently,  as  it  can  be  seen  in  Figure  3.5,  where  the  high  frequency  oscillations  of  vl  are 
linked  to  vortex  shedding.  In  the  spectrum  of  Figure  3.6,  the  high  frequency  is  /  =  1.81,  which 
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reflects  the  frequency  of  continuous  vortex  shedding.  The  lower  /  =  1.28  is  thought  to  reflect  the 
intermittent  feature  of  vortex  shedding.  The  spectrum  of  pressure  fluctuations  exhibit  almost  the 
same  peak  frequency,  e.g.  /  =  1.28  and  /  =  1.79. 


Figure  3.5:  Instantaneous  streamwise  velocity  at  Figure  3.6:  Power  spectrum  density  of 
location  point  PI  ul  at  location  point  PI 


Figure  3.7:  Instantaneous  pressure  at  location  Figure  3.8:  Power  spectrum  density  of 
point  PI  p'  at  location  point  PI 

Figure  3.9  -  3.12  show  the  time  series  of  streamwise  velocity,  pressure  and  their  spectrums  at 
point  P3.  At  this  point,  the  signals  of  streamwise  velocity  and  pressure  are  similar  to  those  at 
point  PI.  The  frequencies  obtained  from  the  velocity  and  pressure  signals  are  almost  the  same. 
The  peak  frequency  is  /  =  1.70,  which  is  close  to  the  vortex  shedding  frequency  at  PI. 

Figure  3.13  -  3:16  show  the  time  series  of  streamwise  velocity  and  pressure  and  their  spectrums 
at  point  P5.  The  instantaneous  velocity  and  pressure  are  more  regular  compared  with  those  at  PI. 
As  a  matter  of  fact,  the  signals  recorded  at  P5  represent  the  reattachment  of  the  primary  vortex, 
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Figure  3.9:  Instantaneous  streamwise  velocity  at  Figure  3.10:  Power  spectrum  density  of 
location  point  P3  u5  at  location  point  P3 


Figure  3.11:  Instantaneous  pressure  at  location  Figure  3.12:  Power  spectrum  density  of 
point  P3  p’  at  location  point  P3 


as  it  is  displayed  in  frames  (f)  (h)  of  Figure  3.2.  The  primary  vortices  shed  from  the  leading  edge 
rotate  clockwise.  When  a  vortex  reattaches  near  point  P5,  which  is  close  to  the  wall,  the  recorded 
streamwise  velocity  decrease,  while  the  vortex  also  causes  the  pressure  to  decrease.  Therefore,  the 
streamwise  velocity  signal  is  in  phase  with  pressure  signal,  as  it  is  shown  in  Figure  3.13  and  Figure 
3.15.  In  Figure  3.13  and  Figure  3.15,  the  low  velocity  and  low  pressure  points  are  corresponding 
to  the  vortex  reattachment.  The  frequencies  obtained  from  the  velocity  and  pressure  signals  are 
almost  the  same.  The  frequency  corresponding  to  the  vortex  reattachment  is  /  =  0.29. 

Point  P9  locates  right  above  the  suction  surface  near  the  trailing  edge.  In  Figure  3.3,  large 
scale  vortex  is  found  near  the  trailing  edge.  The  signals  recorded  at  point  P9  is  used  to  analyze  the 
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Figure  3.13:  Instantaneous  streamwise  velocity  at  Figure  3.14:  Power  spectrum  density  of 
location  point  P5  u’  at  location  point  P5 


Figure  3.15:  Instantaneous  pressure  at  location  Figure  3.16:  Power  spectrum  density  of 
point  P5  p’  at  location  point  P5 


frequency  feature  of  the  trailing  edge  vortex.  The  instantaneous  fluctuations  of  streamwise  velocity 
and  pressure  at  P9  and  their  power  spectrum  are  displayed  in  Figure  3.17  -  Figure  3.20.  Both  the 
velocity  and  pressure  signals  at  point  P9  are  dominated  by  low-frequency  oscillations.  It  is  also 
interesting  to  notice  that  at  this  point,  the  streamwise  velocity  is  opposite  in  phase  to  pressure; 
i.e.,  the  streamwise  velocity  decreases  as  pressure  increases.  A  brief  explanation  is  given  here  based 
on  the  fact  that  the  trailing  vortex  rotates  counterclockwise.  As  the  trailing  vortex  appears,  at 
point  P9  close  to  wall,  the  induced  motion  of  the  vortex  increases  the  streamwise  velocity,  while 
the  pressure  decreases.  At  this  point,  the  peak  frequency  for  velocity  is  /  =  0.24  and  /  =  0.28  for 
pressure. 

Now  we  focus  on  the  next  location  point  Pll,  which  locates  in  the  wake  but  still  close  to  the 
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Figure  3.17:  Instantaneous  streamwise  velocity  at  Figure  3.18:  Power  spectrum  density  of 
location  point  P9  u’  at  location  point  P9 


Figure  3.19:  Instantaneous  pressure  at  location  Figure  3.20:  Power  spectrum  density  of 
point  P9  p’  at  location  point  P9 

trailing  edge.  High-frequency  oscillations  appear  on  streamwise  velocity  and  pressure  displayed  in 
Figure  3.21  and  3.23,  respectively.  The  high-frequency  oscillations  associated  with  the  small-scale 
vortex  structures  appear  from  the  trailing  edge  along  the  edge  the  large-scale  trailing  edge  vortex, 
as  shown  in  frames  (i),  (j)  (k)  of  Figure  3.3.  The  high-frequency  parts  are  modulated  by  the 
low-frequency  signals  corresponding  to  the  large-scale  trailing  edge  vortex. 

Point  P12  locates  at  downstream  of  Pll,  the  signals  recorded  at  initial  stage  at  P12  correspond¬ 
ing  to  frame  (a)  to  (c)  in  Figure  3.3  are  displayed  from  Figure  3.25  to  Figure  3.28.  The  regular 
pattern  of  the  velocity  and  pressure  indicate  the  typical  Kelvin-Helmholtz  type  instability. 
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Figure  3.21:  Instantaneous  streamwise  velocity  at  Figure  3.22:  Power  spectrum  density  of 
location  point  Pll  u’  at  location  point  Pll 


Figure  3.23:  Instantaneous  pressure  at  location  Figure  3.24:  Power  spectrum  density  of 
point  Pll  p’  at  location  point  Pll 

3.5  Conclusions 

Direct  numerical  simulation  is  carried  out  by  solving  the  full  Navier-Stokes  equations  in  generalized 
curvilinear  coordinates  to  study  the  separation  flow  around  a  NACA  0012  airfoil  at  large  angle  of 
attack.  By  using  a  fourth-order  centered  compact  scheme  for  spatial  discretization,  the  small-scale 
vortical  structures  are  resolved,  which  will  dissipate  if  low-order  numerical  schemes  are  used.  Non¬ 
reflecting  boundary  conditions  are  imposed  at  the  far-field  and  outlet  boundaries  to  avoid  possible 
non-physical  wave  Reflection. 

The  DNS  results  clearly  describe  the  flow  separation  process  at  the  upper  surface  of  the  airfoil. 
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Figure  3.25:  Early  stage  of  the  instantaneous 
streamwise  velocity  at  location  point  P12 
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Figure  3.26:  Power  spectrum  density  of 
signal  u’  in  Figure  3.25 
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Figure  3.27:  Early  stage  of  the  instantaneous  pres-  Figure  3.28:  Power  spectrum  density  of 
sure  at  location  point  P12  signal  p’  in  Figure  3.27 


The  phenomena  of  the  leading-edge  separation,  vortex  shedding,  vortex  merging,  vortex  pairing, 
and  formation  and  shedding  of  large-scale  trailing  edge  vortex  axe  displayed  and  discussed  in  detail. 
The  small-scale  vortices  associated  with  the  Kelvin-Helmholtz  instability  are  also  observed  along  the 
shear  layer  near  the  trailing-edge.  These  phenomena  are  in  good  agreement  with  the  experimental 
results  obtained  by  Shih  et  al.  (1992,  1995). 
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Chapter  4 

Numerical  Simulation  of  Flow  Instability  around  a 
Delta  Wing 


4.1  Introduction 

Recent  developments  in  aerospace  technology  have  revived  the  interest  to  the  study  of  flow  separa¬ 
tions  around  an  aircraft  maneuvering  dynamic  operations.  Understanding  of  the  complex  separated 
vortical  flow  is  crucial  to  the  aerodynamic  design  of  modern  aircraft.  Vortical  structures,  which 
develop  over  the  leading-edge  extension,  slender  fore-body,  and  main  wing,  may  have  severe  effect 
on  the  aerodynamic  characteristics  and  performance  of  modern  fighter  aircraft. 

A  flat-plate  delta  wing  with  sharp  leading-edge  provides  a  simple  configuration  to  investigate 
the  development  of  the  vortical  structures.  Both  experimental  and  computational  results  have 
shown  that  the  flow  over  the  suction  side  of  a  delta  wing  at  a  fixed  angle  of  attack  is  dominated 
by  a  pair  of  counter-rotating  vortices,  i.e.  the  leading-edge  primary  vortices.  These  vortices  are 
formed  as  a  result  of  the  rolling-up  of  the  vortex  sheet  shedding  from  the  leading-edge.  The  flow 
induced  by  the  leading-edge  vortices  separates  near  the  wing  surface  and  iorms  a  pair  of  oppositely 
rotating  secondary  vortices.  The  size  and  strength  of  the  leading-edge  vortices  increase  with  the 
angle  of  incidence,  resulting  in  a  substantial  nonlinear  lift  increment.  But  the  maximum  lift  of  a 
delta  wing  is  limited  by  a  phenomenon  known  as  vortex  breakdown (Visser  &  Nelson  1993). 

In  the  experimental  study  using  dye  visualization  carried  out  by  Gad-el-Hak  and  Balckwelder 
(1985),  small-scale  vortices  were  observed  shedding  from  the  leading-edge  and  feeding  into  the 
rolling-up  process  associated  with  the  large-scale  leading-edge  vortices.  In  their  study,  the  small- 
scale  vortices  were  paired,  and  were  believed  as  the  origination  of  the  classical  leading-edge  vortices. 
Payne  et  al.  (1988)  used  smoke  flow  visualization  and  laser  sheet  technique  to  study  the  vortical 
flow  field  above  the  delta  wing  at  high  angles  of  attack.  Two  types  of  vortex  breakdown  were 
testified,  i.e.  the  bubble  mode  and  the  spiral  mode.  In  the  same  experiment,  static  small-scale 
vortical-like  structures  were  found  in  the  shear  layer  of  a  delta  wing  with  a  85°  sweep  angle  and 
a  40°  angle  of  attack.  The  growth  of  these  structures  was  found  to  be  similar  to  the  evolution  of 
the  classic  Kelvin- Helmholtz  instability.  In  this  experiment,  the  pairing  of  the  small-scale  vortices 
was  not  observed.  .The  recent  experimental  work  of  Rieley  &  Lowson  (1998)  revealed,  using  flow 
visualization  and  hot-wire  measurement,  the  existence  of  static  small  vortical  structures  in  the 
free  shear  layer  shedding  from  the  sharp  leading-edge  of  a  delta  wing.  A  local  three-dimensional 
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Kelvin-Helmholtz-type  instability  was  suggested  by  the  authors  for  the  formation  of  these  vortical 
structures  in  the  free  shear  layer.  Similar  vortical  structures  were  also  observed  in  the  investigations 
of  Cipolla  &  Rockwell  (1998),  where  small-scale  concentrations  of  vorticity  form  near  the  leading- 
edge  of  a  rolling  delta  wing.  These  vortices  appear  to  evolve  in  a  coupled  fashion,  which  has  been 
considered  as  the  wake-like  instability. 

Numerical  simulations  of  vortex  breakdown  above  a  stationary  sharp  edged  delta  wing  over 
a  range  of  angles  of  attack  were  carried  out  by  Modiano  &  Murman  (1994).  Their  computation 
was  based  on  an  Euler  solver  with  adaptive  mesh  system.  The  spiral  form  of  vortex  breakdown 
was  observed  without  the  emergence  of  the  small-scale  vortical  structures  inside  the  shear  layer. 
In  the  numerical  investigation  by  Argwal  et  al  (1992),  the  well-known  Euler/Navier-Stokes  code 
CFL3D  was  used  to  simulate  the  leading-edge  vortex  breakdown  of  a  low-speed  flow  on  a  flat-plate 
delta  wing  with  sharp  leading-edges.  Although  the  vortex  breakdown  positions  obtained  from  the 
computation  were  reported  in  good  agreement  with  experimented  data,  the  small-scale  vortices  were 
not  observed,  which  could  be  attributed  to  the  lack  of  numerical  resolution/accuracy.  A  numerical 
investigation  of  the  unsteady  vortex  structure  over  a  76°  sweep  wing  at  20.5°  angle  of  attack  was 
carried  out  by  Gordnier  &  Visbal  (1994).  Their  numerical  calculation  indicated  that  the  small-scale 
vortical  structures  emanated  from  the  leading-edge  was  brought  on  by  the  Kelvin-Helmholtz-type 
instability.  Pairing  of  the  small  vortices  was  not  observed  in  the  computational  results. 

The  intention  of  present  work  is  to  study  the  mechanism  of  vortex  breakdown  above  a  slender 
flat-plate  delta  wing  with  sharp  leading-edges  at  a  fixed  angle  of  attack.  Direct  numerical  simulation 
is  employed  to  give  a  detailed  description  of  flow  instability  and  vortex  shedding  near  the  leading- 
edge  of  the  delta  wing. 


4.2  Governing  Equations 


The  three-dimensional  compressible  Navier-Stokes  equations  in  generalized  curvilinear  coordinates 
(£,  rj,Q  are  written  in  conservative  forms: 


ldQ  djE-Ey)  djF-FV)  djG-Gy)  . 

J  dt  d£  &rj 


The  flux  vectors  for  compressible  flow  are 
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where;  J  is  Jacobian  of  the  coordinate  transformation,  and  ,  £y,  ,  Vx ,  rly ,  J?* i  Cx  t  Cy  >  Cz  are  coordi¬ 

nate  transformation  metrics,  r^’s  and  g^’s  are  the  viscous  stress  and  the  heat  flux,  respectively. 

In  Eq.  (4.1),  second  order  Euler  Backward  scheme  is  used  for  time  derivatives,  and  the  fully 
implicit  form  of  the  discretized  equations  is: 

3Q"+i  _  4 Q"  +  <7»-i  9(£"+1  -  E%+1) 


Gv  =  ~ F 
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dd 
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=  0 


gn+1  =  qp  +  jqp 


(4.2) 

(4.3) 


Qn+1  is  estimated  iteratively  as: 
where, 

5QP  =  Qp+1  -  QP  (4.4) 

At  step  p  =  0,  Qv  =  Qn]  as  SQP  is  driven  to  zero,  Qv  approaches  Qn+l  ■  Flux  vectors  are  linearized 
as  follows: 


En+1  «  EP  +  AP8QP 
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_F"+1  «  fP  +  BP5QP 
Gn+l  ^Gp  +  CpSQp 


So  that  Eq.  (4.2)  can  be  written  as: 


gj  +  A  tJ{D^A  +  D„B  +  D$C)]5QP  =  R 

£ 


where  R  is  the  residual: 


R  =  -^Qp-2Qn  +  ^Qn~1)-AtJ[(D^E-Ev) 

+  DV{F-Fv)+DAG-Gv)]p 

D^Dj^D^  represent  partial  differential  operators,  and  A,  B,  C  are  the  Jacobian  matrices  of  flux 


vectors: 


dE  dF  dG 

A  dQ ’  B  dQ'  dQ 


dQ'  dQ'  dQ 

The  right  hand  side  of  Eq.  (4.6)  is  discretized  using  fourth-order  compact  scheme  (Lele,  1992) 
for  spatial  derivatives,  and  the  left  hand  side  of  the  equation  is  discretized  following  LU-SGS  method 
(Yoon  &  Kwak,  1992).  In  this  method,  the  Jacobian  matrices  of  flux  vectors  are  split  as: 

A  =  A+  +  A“,  B  =  B+  +  B~,  C  =  C+  +  C~ 


where, 


A±  =  ~[A  ±  rAI\ 
B±  =  \[B±tbI] 
C±  =  \[C±tcI\ 


ta  —  «raaa;[|A(A)|]  +  v 
tb  =  /cma:r[|A(i?)|]  +  v 
rc  =  nmax[\X(C)\]  +  v 

where  A(A),  A(£),  A(C)  are  eigenvalues  of  A,  B ,  C  respectively,  k  is  a  constant  greater  than  1.  v  is 
taken  into  account  for  the  effects  of  viscous  terms,  and  the  following  expression  is  used: 

V  =  Wax[(7-l)M^ePr’  3  7fJ 

The  first-order  upwind  finite  difference  scheme  is  used  for  the  split  flux  terms  in  the  left  hand 
side  of  Eq.  (4.6).  This  does  not  effect  the  accuracy  of  the  scheme.  As  the  left  hand  side  is  driven 
to  zero,  the  discretization  error  will  also  be  driven  to  zero.  The  finite  difference  representation  of 
Eq.  (4.6)  can  be  written  as: 

[|/  +  AtJ(rA  +  rB  +  rc)I}6QpiJ>k  =  R?iJik 
-AtJ  [  ASQ?+1Jik-A+6Q?_1Jik 
+  BS^+lik-B+5^_ltk 

+  C-6<yij,k+i-C+6QPij,k-i  1  (4-9) 
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In  LU-SGS  scheme,  Eq.  (4.9)  is  solved  by  three  steps.  First  initialize  SQ°  using 

SQl,k  =  [\l  +  A  tJ(rA  +rB  +  rc)I]~lRl,k 

In  the  second  step,  the  following  relation  is  used: 

SQid,k  =  6Qij,k  +  [p  +  AtJ^  +  rB  +  re)/]-1 

x[A  tJ(A+6QUj,k  +  B+SQ;j_hk  +  C+JQ^-i)] 

For  the  last  step,  5QP  is  obtained  by 

S&ijj,  =  6Q'ij,k  -  [p  +  MJ(rA  +  TB  +  r c)]~l 
x[A tJ(A-6<?i+1>j>k  +  B-8Q?J+l>k  +  C-5Cyi)jM1)] 

The  sweeping  is  performed  along  the  planes  of  i+j +k  —  const ,  i.e.  in  the  second  step,  sweeping 
is  from  the  low-left  corner  of  the  grid  to  the  high-right  corner,  and  then  vice  versa  in  the  third  step. 

In  order  to  depress  numerical  oscillation  caused  by  central  difference  scheme,  spatial  filtering  is 
used  instead  of  artificial  dissipation.  Implicit  sixth-order  compact  scheme  for  space  filtering  (Lele, 
1992)  is  applied  for  primitive  variables  u,  v,  w ,  p,p  after  each  time  step. 

For  subsonic  flow,  u,  v,w,  T  are  prescribed  at  the  inflow  boundary,  p  is  obtained  by  solving 
modified  N-S  equation  based  on  the  characteristic  analysis.  On  the  far  field  and  out  flow  boundary, 
non- reflecting  boundary  conditions  are  applied.  Adiabatic,  non-slipping  conditions  are  used  for  the 
wall  boundary.  All  equations  of  boundary  conditions  are  solved  implicitly  with  internal  points. 
Specific  details  of  boundary  treatment  can  be  found  in  Jiang  et  al  (1999). 


4.3  Computational  Details 

Direct  numerical  simulation  has  been  implemented  to  investigate  the  compressible  flow  separation 
around  a  slender  delta  wing.  The  geometry  of  the  delta  wing,  taken  from  the  experimental  work 
of  Rieley  &  Lowson  (1998),  is  shown  in  Figure  4.1.  The  sweep  angle  denoted  by  A  is  80°  and  the 
leading-edge  angle  denoted  by  o  is  30°.  The  chord  length  is  taken  as  the  characteristic  length  L, 
such  that  the  non-dimensional  chord  length  is  c  =  I.OjL.  The  non-dimensional  thickness  of  the  delta 
wing  is  h  =  0.024 L.  The  freestream  velocity  t/oo  is  the  characteristic  velocity. 

An  H-C  type  mesh  system  for  a  half-plane  model  of  the  delta  wing  is  used  based  on  the  as¬ 
sumption  that  the  flow  is  symmetrical  to  the  the  half-plane.  The  mesh  is  H-type  in  the  meridian 
section  and  C-type  in  the  cross  section.  An  elliptic  grid  generation  method,  first  proposed  by 
Spekreijse  (1995),  is  used  to  generate  the  three-dimensional  grids.  This  method  is  based  on  a 
composite  mapping,  which  is  consisted  of  a  nonlinear  transfinite  algebraic  transformation  and  an 
elliptic  transformation.  The  grids  are  orthogonal  on  the  delta  wing  surface.  The  sharp  leading-edge 
is  approximated  by  a  round  edge  with  a  small  radius  of  1.0  x  10 ”3L,  while  in  the  experiment  of 
Rieley  &  Lowson  (1998),  the  average  thickness  of  the  leading-edge  was  0.12  mm,  which  was  approx¬ 
imately  2.55  x  10-4L.  Computations  are  carried  out  on  two  grid  systems,  i.e.  the  low  resolution 
mesh  with  140  x  70  x  70  grid  nodes  and  the  high  resolution  mesh  with  180  x  150  x  70  grid  nodes, 
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Figure  4.1:  Schematic  of  the  delta  wing 


Figure  4.2:  H-C  type  grid  around  a  85°  sweep  deH;.  wing 

where  the  sequence  of  numbers  is  corresponding  to  the  axial,  the  spahwise  and  the  wall-normal 
direction,  respectively.  An  example  of  the  three-dimensional  grid  is  displayed  in  Figure  4.2. 

The  parallel  version  of  the  DNS  code  based  on  the  Message  Passing  Interface  (MPI)  has  been 
developed  to  accelerate  the  computation.  Although  massive  data  exchanges  axe  required  for  com¬ 
puting  the  derivatives  with  the  compact  finite  difference  scheme,  the  speedup  is  still  substantial. 
On  a  SGI  Origin  2000  computer,  the  performance  of  the  MPI  code  is  superior  to  the  serial  code, 
which  is  compiled  using  the  automatic  parallelization  option  provided  by  Fortran  90  compiler.  The 
comparison  are  displayed  in  Figure  4.3,  where  the  speedup  parameter  S(n,p)  is  defined  as  the  ratio 
of  the  runtime  of  a  serial  solution  to  a  problem  to  the  parallel  runtime.  Linear  speedup  has  been 
achieved  for  MPI  code  running  on  4,  6,  and  15  processors.  The  final  parallel  computation  for  the 
higher  Reynolds  number  case  with  a  mesh  of  180  x  150  x  70  is  conducted  using  10  processors  on  a 
SGI  Origin  2000  computer. 
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Figure  4.3:  Speedup  S{n,p)  of  the  MPI  code  compared  with  the  serial  code  complied  with  automatic 
parallelization  option 

4.4  Results  and  Discussions 

4.4.1  Low  Reynolds  number  case 

All  results  presented  here  are  obtained  from  direct  numerical  simulation  of  flow  around  a  85° 
sweep  delta  wing  with  a  flat-plate  suction  surface,  which  has  been  introduced  in  Section  3.  The 
angle  of  attack  is  fixed  at  12.5°.  The  free-stream  Mach  number  is  Ma  =  0.1.  For  the  lower 
resolution  simulation,  the  Reynolds  number  based  on  the  chord  length  and  the  free-stream  velocity 
is  Rec  =  5  x  104.  No  initial  and  boundary  disturbance  axe  imposed  in  the  simulation. 

The  contours  of  the  axial  vorticity  on  selected  cross  sections  are  displayed  in  Figure  4.4.  It 
is  quite  cleax  that  a  pair  of  counter-rotating  vortices,  so  called  the  leading-edge  primary  vortices, 
appears  over  the  suction  side  of  the  delta  wing.  These  vortices  form  as  a  result  of  flow  separation 
and  the  rolling-up  of  the  vortex  sheet  shedding  from  the  leading-edge.  The  primary  vortices  are 
steady  and  stable  in  this  low  Reynolds  number  case.  The  primary  vortices  are  composed  of  a 
pair  of  counter-rotating  oblique  vortex  tubes  starting  from  the  apex  of  delta  wing,  from  a  three- 
dimensional  point  of  view.  Beneath  the  primary  vortices,  near  the  upper  surface  of  the  delta  wing, 
the  secondary  vortices,  which  have  an  opposite  rotating  direction  to  the  primary  vortices,  are  formed 
as  a  result  of  the  spanwise  outflow  induced  by  the  primary  vortex.  Figure  4.5  shows  the  three- 
dimensional  instantaneous  streamlines  starting  from  vicinity  of  the  leading-edge.  The  streamlines 
also  reveals  the  existence  of  the  cone-shaped  primary  vortices.  The  computational  results  are  in 
good  agreement  with  the  experimental  results  of  Riley  &  Lowson  (1998).  During  the  computation 
for  the  high  Reynolds  number  case,  the  flow  becomes  unsteady  and  small-scale  vortical  structures 
keep  shedding  from  the  leading-edge.  In  the  experiment  of  Riley  &  Lowson  (1998),  flow  instability 
was  observed  when  the  Reynolds  number  was  raised  above  Rec  =  100, 000.  In  order  to  study  the 
flow  instability  near  the  leading  edge,  the  numerical  simulation  with  a  higher  resolution  and  higer 
Reynolds  number  has  been  implemented,  and  the  results  are  discussed  in  next  section. 
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Figure  4.4:  Contours  of  the  axial  vorticity 
on  selected  cross  sections,  angle  of  attack 
a  =  12.5°,  Re  =  5  x  104,  Mo  =  0.1 


Figure  4.5:  Three-dimensional  streamline 
above  the  suction  surface  of  a  85°  sweep 
delta  wing  at  an  angle  of  attack  a  =  12.5°, 
Re  =  5  x  104,  Ma  =  0.1 


4.4.2  High  Reynolds  number  case 

At  higher  Reynolds  number,  i.e.  Rec  =  1.96  x  105,  which  is  the  same  as  the  one  chosen  by  Riley  &; 
Lowson  (1998)  in  their  experimental  study,  where  the  flow  instability  occurs  near  the  leading  edge 
of  the  delta  wing.  In  order  to  capture  the  small  vortical  structures  observed  in  the  experiment,  the 
numerical  simulation  is  accomplished  on  a  mesh  with  a  higher  resolution  of  180  x  150  x  70.  During 
the  simulation,  flow  instability  and  periodic  shedding  of  small  vortical  structures  are  observed. 
Since  there  is  no  disturbance  imposed  as  the  initial  or  boundary  condition  for  the  computation,  the 
unstable  behavior  presented  by  the  flow  in  the  computational  results  are  rather  intrinsic. 

The  distributions  of  the  instantaneous  axial  vorticity  on  various  cross  sections  are  shown  in 
Figure  4.6.  Compared  with  the  low  Reynolds  number  results  of  Figure  4.4,  the  flow  is  still  dominated 
by  a  pair  of  primary  vortices.  But  instability  appears  at  the  leading-edge  of  delta  wing,  where  vortex 
shedding  is  observed.  On  the  suction  surface  near  the  leading-edge,  the  secondary  vortices  are  still 
visible  in  this  figure. 

In  Figure  4.7  the  contours  of  axial  vorticity  at  different  time  on  a  cross  section  at  x  =  0.87L 
are  displayed  through  (a)  to  (h),  each  frame  is  corresponding  to  a  snapshot  of  a  two-dimensional 
flow  field  at  a  certain  time.  Flow  instability  is  quite  obvious  in  these  figures.  The  primary  vortex 
deforms  compared  to  the  low  Reynolds  number  case.  The  flow  pattern  inside  the  primary  vortex 
resembles  that  of  the  spiral  instability  mode,  which  presented  occasionally  in  the  experiment  of 
Rieley  &  Lowson(1998).  Two  strong  shear  layers  are  visible  through  the  concentration  of  axial 
vorticity  contours.  The  first  one  is  the  leading-edge  shear  layer  whose  axial  vorticity  is  positive 
(shown  in  light  color  in  Figure  4.7),  which  wraps  the  leading-edge  corner  from  below  and  feeds 
into  the  primary  vortex.  The  other  one  lies  between  the  primary  vortex  and  the  upper  surface 
of  the  delta  wing  and  has  a  negative  axial  vorticity  (shown  in  dark  color  in  Figure  4.7),  which  is 
associated  with  the  secondary  vortex.  Therefore,  the  shear  layer  below  the  primary  vortex  is  also 
called  the  secondary  shear  layer.  As  it  will  be  discussed  later,  both  the  leading-edge  shear  layer 
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Figure  4.6:  The  instantaneous  distributions  of  the  axial  vorticity  on  various  cross  sections.  Angle 
of  attack  a  =  12.5°,  Re  =  1.96  x  105,  Ma  —  0.1 


and  the  secondary  shear  layers  are  related  to  the  instability  and  vortex  shedding  procedure  near 
the  leading-edge. 

Among  the  small-scale  vortical  structures  shedding  from  the  leading-edge,  there  are  two  types  of 
vortices,  distinguished  by  the  direction  of  rotation  or  by  the  sign  of  axial  vorticity.  Those  vortices 
whose  rotating  direction  is  the  same  as  the  primary  vortex  are  named  as  the  A-family  vortices, 
which  are  corresponding  to  a  positive  axial  vorticity  component.  The  vortices  rotating  in  the 
opposite  direction  as  the  primary  vortex  are  called  the  B- family  vortices  and  have  a  negative  axial 
vorticity  component.  The  A-family  vortices  are  much  stronger  than  the  B-family  vortices,  which 
can  be  recognized  from  the  contours  of  the  axial  vorticity  in  Figure  4.7. 

In  Figure  4.7(a),  a  bulge  is  observed  on  the  leading-edge  shear  layer.  The  bulge  grows  in  size 
as  it  moves  away  from  the  leading  edge,  as  shown  in  Figure  4.7(b),  (c),  and  (d).  This  process  is 
corresponding  to  the  generation  and  shedding  of  the  A-family  vortex.  Because  the  B-family  vortices 
are  very  weak,  the  ■  bedding  process  of  B-family  vortices  is  not  clear  in  Figure  4.7.  A  more  detailed 
study  reveals  tha*  the  B-family  vortex  comes  from  the  shedding  of  the  secondary  vortex. 


(g)  t  =  56.37L/U, o, 


(h)  t  =  56.51L/{7a 


Figure  4.7:  Contoiirs  of  axial  vorticity  of  different  time  on  a  cross  section  at  x  =  0.87L.  Angle  of 
attack  a  =  12.5°,  Re  —  1.96  x  105,  Ma  =  0.1 
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In  Figure  4.7,  frame  (a)  resembles  frame  (g).  So  frames  (a)  through  (g)  show  one  period 
of  vortex  shedding.  The  small-scale  vortical  structures  shedding  from  the  leading-edge  are  in  turn 
captured  by  the  primary  vortex  and  feeding  into  the  rolling-up  process  of  the  primary  vortex.  After 
leaving  the  leading-edge,  both  A-  and  B- family  vortices  experience  severe  deformation  as  they  are 
stretched  and  captured  by  the  primary  vortex.  The  flow  pattern  inside  the  core  of  the  primary 
vortex  resembles  that  of  the  spiral  instability  mode,  which  is  also  observed  in  the  experiment  (Rieley 
&  Lowson,  1998).  In  the  numerical  simulation,  these  small-scale  vortical  structures  dissipate  quickly 
as  they  are  leaving  the  leading-edge  and  entering  the  central  region  of  the  primary  vortex.  But  in 
the  experiment  the  spiral  instability  can  evolve  into  transition.  The  fast  dissipation  of  the  spiral 
mode  in  the  numerical  simulation  can  be  attributed  to  the  insufficiency  of  resolution.  Because  grids 
are  clustered  near  the  wall  and  near  the  leading-edge  of  delta  wing,  in  the  areas  far  away  from  the 
wall  the  resolution  is  relatively  low.  The  future  effort  will  be  devoted  to  increasing  the  resolution 
where  the  primary  vortex  locates. 

Because  the  small-scale  vortical  structures  are  shedding  continuously  from  the  leading-edge, 
the  area  near  the  leading-edge  is  of  particular  interesting  in  present  work.  The  detailed  pictures 
of  vortex-shedding  near  the  leading-edge  is  shown  in  Figure  4.8,  where  the  limiting  streamline  and 
contours  of  axial  vorticity  of  various  instance  on  a  cross-section  at  x  =  0.87 L  axe  displayed.  Through 
(a)  to  (h)  in  Figure  4.8,  the  pattern  of  limiting  streamline  exhibits  a  periodic  feature.  Actually, 
figures  (a)  to  (g)  fit  in  one  period  of  variation.  In  Figure  4.8(a),  there  is  a  secondary  vortex 
attaching  on  the  upper  surface  of  delta  wing.  The  strong  leading-edge  shear  layer  is  shown  by  dark 
color  of  the  contours  of  axial  vorticity.  Near  the  leading-edge,  the  shear  layer  is  concentrated  in  a 
narrow  area.  In  Figure  4.8(b),  at  t  =  55.68L/f7oo>  a  small  vortex  shown  by  the  limiting  streamline 
appears  near  the  leading-edge  over  the  free  shear  layer.  In  the  same  picture,  a  small  bulge  appears 
on  the  shear  layer.  The  generation  of  this  small  vortex  can  be  attributed  to  the  Kelvin-Helmholtz 
instability.  Therefore,  the  small  vortex  is  named  as  the  Kelvin-Helmholtz  (K-H)  type  vortex.  At 
the  same  time,  the  secondary  vortex,  which  was  attaching  on  the  wing  surface  at  t  =  55.54L/J7QO, 
moves  away  from  the  wall.  As  the  K-H  type  vortex  grows,  the  secondary  vortex  is  pushed  further 
away  from  the  wall.  From  t  =  55.8lL/Uoo  to  t  =  55.95L/f/oo,  (Figure  4  °(o)  to  (d))  the  secondary 
vortex  moves  upward  and  begins  to  separate  from  the  wall,  which  is  corresponding  to  the  B-family 
vortex,  whose  rotating  direction  is  opposite  to  the  primary  vortex.  Therefore,  the  B-family  vortex 
comes  from  the  shedding  of  the  secondary  vortex.  The  generation  of  the  leading-edge  K-H  type 
vortex  also  causes  the  deformation  of  the  shear  layer,  which  is  visible  from  the  contours  of  the  axial 
vorticity  in  Figure  4.8(b),  (c),  and  (d).  The  bulge  on  the  contours  of  axial  vorticity  is  corresponding 
to  the  K-H  type  vortex.  In  Figure  4.8(d),  the  secondary  vortex  almost  disappears  and  the  K-H  type 
vortex  is  still  attached  to  the  leading-edge.  In  Figure  4.8(e)  and  (f),  the  K-H  type  vortex  grows  in 
size  until  it  reaches  the  edge  of  the  primary  vortex.  Another  vortex  appears  at  the  same  location  of 
the  secondary  vortex,  actually  it  is  a  new  secondary  vortex.  The  K-H  type  vortex  gradually  moves 
upward  and  sheds  from  the  leading-edge,  and  comes  out  to  be  the  A- family  vortex,  whose  rotating 
direction  is  the  same  as  the  primary  vortex.  It  is  obviously  that  the  A-family  vortex  originates 
from  the  K-H  type  leading-edge  vortex.  The  last  two  frames  are  the  periodic  repeating  of  frames 
(a)  and  (b)  in  Figure  4.8. 

The  vortex-shedding  near  the  leading-edge  is  a  periodic  process.  The  interaction  between  the 
secondary  vortex  and  the  leading-edge  shear  layer  generates  a  K-H  type  vortex.  As  this  K-H  type 
vortex  grows,  the  induced  flow  pushes  the  secondary  vortex  away  form  the  wall,  and  ultimately 
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leads  to  the  shedding  of  the  B-family  vortex.  The  K-H  type  vortex  grows  in  size  as  the  secondary 
vortex  shows  up  again  near  the  wall.  The  induced  flow  pushes  the  K-H  type  vortex  away  form  the 
wall  and  leads  to  the  shedding  of  the  A-family  vortex.  So  the  A-family  vortex  originates  from  the 
Kelvin-Helmholtz  instability  of  shear  flow  near  the  leading-edge.  The  B-family  vortex  originates 
from  the  secondary  vortex.  The  period  of  vortex  shedding  is  between  O.SQL/Uoo  and  0.98Z//f7oo- 
The  scale  of  the  K-H  type  leading-edge  vortex  and  the  secondary  vortex  is  about  0.005L. 

The  interpretation  of  the  above  phenomena  is  based  on  the  Kelvin-Helmholtz  instability  6f 
cross-sectional  two-dimensional  flow.  Considering  many  cross-sections  simultaneously,  the  period 
of  vortex  shedding  is  the  same,  there  is  only  phase  difference  between  one  cross-section  and  the 
other.  From  a  three-dimensional  point  of  view,  the  A-  and  B-family  vortices  become  vortex  tubes, 
which  are  oblique  to  the  axial  direction. 
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(g)  t  =  56.37L/E/oo  (h)  t  =  56.51L/J7oo 

Figure  4,8:  Limiting  streamline  and  contours  of  axial  vorticity  of  different  time  on  a  cross  section 
at  x  =  0.87 L.  Angle  of  attack  a  =  12.5°,  Re  =  1.96  x  105,  Ma  =  0.1 


The  time  series  of  three  components  of  the  instantaneous  velocity  at  a  location  near  the  leading- 
edge  {x  =  0.87L,  y  =  0.076L,  z  =  0.00942,)  are  recorded  and  shown  in  Figure  4.9,  4.10,  and  4.11. 
This  probe  point  locates  on  the  cross-section  shown  in  Figure  4.7  and  Figure  4.8,  so  that  the  velocity 
signals  can  be  interpreted  in  accordance  with  the  two-dimensional  vortex  shedding  pictures.  The 
signals  of  the  three  components  of  velocity  are  all  periodic  functions  of  time.  The  axial  velocity 
u  has  two  local  maximums  and  two  local  minimunis  within  a  period.  There  are  only  one  local 
maximum  and  one  local  minimum  within  a  period  for  the  signals  of  the  spanwise  velocity  v  and 
the  vertical  velocity  w.  The  phase  difference  between  v  and  w  is  approximately  7r/2,  which  can  be 
interpreted  as  a  result  of  the  small-scale  vortical  structure  passing  through  the  probe. 


Figure  4.9:  Instantaneous  axial  velocity  at  a  location  of  x  =  0.87 2,,  y  =  0.0762,,  z  —  0.00942, 


Figure  4.10:  Instantaneous  spanwise  velocity  at  a  location  of  x  —  0.87L,  y  =  0.0761/,  z  —  0.00942, 


Figure  4.11:  Instantaneous  vertical  velocity  at  a  location  of  x  =  0.872,,  y  =  0.0762,,  z  =  0.00942, 

In  order  to  compare  the  axial  velocity  signal  with  the  vortex  and  shedding  pictures  in  Figure  4.8, 
a  part  of  Figure  4.9  has  been  enlarged  and  shown  in  Figure  4.12,  where  those  points  with  the  same 
time  as  the  frames  of  Figure  4.8  have  been  marked  with  the  same  letter  through  (a)  to  (g).  In 
Figure  4.12,  one  period  starts  at  point  (a)  and  ends  at  point  (g).  Compared  with  Figure  4.8,  it  has 
been  found  that  the  local  minimum  at  point  (d)  with  a  smaller  axial  velocity  value  is  corresponding 
to  the  B-family  vortex.  As  the  B-family  vortex  moves  through  the  probe  point,  as  shown  in 
Figure  4.8(c)  and  (d),  the  axial  velocity  u  decreases  from  point  (c)  to  point  (d)  and  reaches  its 
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local  minimum  in  Figure  4.12.  Then  the  axial  velocity  recovers  as  the  B-family  vortex  leaves  the 
probe  point.  It  is  followed  by  the  shedding  of  an  A-family  vortex  from  the  leading-edge.  Before  the 
central  part  of  the  A-family  vortex  reaches  the  probe  point,  the  axial  velocity  u  signals  recorded 
by  the  probe  increases  from  point  (e)  to  (f)  in  Figure  4.12.  Then  it  decreases  again  from  (f)  to 
(g)  as  the  core  of  the  A-family  vortex  moves  tlirough  the  probe.  Therefore,  the  local  minimum  at 
point  (g)  is  associated  with  the  A-family  vortex.  The  local  minimum  at  point  (g)  has  a  relatively 
larger  value  of  axial  velocity  compared  with  the  local  minimum  at  point  (d),  which  is  related  to 
the  B-family  vortex.  The  center  of  both  A-  and  B-family  vortices  is  low-momentum  region.  Since 
the  B-family  vortex  originates  from  the  shedding  of  the  secondary  vortex  near  the  upper  surface 
of  the  delta  wing,  and  it  brings  fluid  with  lower  axial  velocity,  the  central  part  of  the  B-family 
vortex  has  a  lower  momentum.  The  A-family  vortex  comes  from  the  shedding  of  leading-edge  K-H 
type  vortex.  It  brings  fluid  from  the  free  shear  layer,  which  has  a  relatively  larger  momentum.  In 
Figure  4.12,  the  local  maximum  at  point  (b)  is  corresponding  to  frame  (b)  of  Figure  4.8,  where  both 
the  A-  and  B-family  vortices  have  not  separated  from  the  delta  wing.  The  signals  of  the  spanwise 
and  the  vertical  velocity  can  be  interpreted  in  a  similar  way.  In  Figure  4.11,  the  local  maximum 
is  corresponding  to  the  passing  of  a  B-family  vortex,  where  the  vertical  velocity  is  positive.  The 
local  minimum  is  corresponding  to  the  passing  of  an  A-family  vortex,  where  the  vertical  velocity  is 
negative.  Thus  the  period  of  velocity  signals  reflects  the  elapsed  time  at  which  the  A-  and  B-family 
vortices  are  shedding  from  the  delta  wing.  Thus  the  period  of  vortex-shedding  can  be  measured 
as  the  distance  between  the  peaks  of  local  maximums  or  local  minimums  on  the  signals  of  three 
velocity  components. 


Figure  4.12:  Instantaneous  axial  velocity  at  location  x  =  0.87L,  y  ==  0.076L,  z  =  0.0094L 

In  order  to  estimate  the  vortex-shedding  period  more  accurately,  power  spectrums  of  velocity 
fluctuations  are  calculated  based  on  the  velocity  signals  and  shown  in  Figure  4.13,  4.14,  and  4.15. 
The  velocity  fluctuation  is  calculated  as  the  difference  between  the  instantaneous  velocity  and 
the  temporal-averaged  mean  velocity.  There  are  two  peaks  in  the  spectrum  of  the  axial  velocity 
fluctuation  shown  in  Figure  4.13.  The  frequency  of  the  first  peak  is  1.086f/oo/£  and  2.31£/oo/£ 
for  the  second  peak.  There  is  only  one  peak  in  the  spectrum  of  the  spanwise  and  vertical  velocity 
fluctuation.  The  peak  frequency  is  1.086[/oo/I',  which  is  the  same  as  the  first  peak  of  the  v! 
spectrum  in  Figure  4.13.  This  peak  frequency  value  /  =  1.086^7oo/J5  represents  the  frequency  of 
vortex  shedding,  the  corresponding  period  is  T  =  0.9208L/Uoo .  The  higher  frequency  in  Figure  4.13, 
corresponding  to  a  time  period  of  T  =  OASSL/Uqq,  reflects  the  elapsed  time  between  the  shedding 
of  a  A-family  vortex  and  a  B-family  vortex,  which  is  approximately  half  the  period  of  the  shedding 
of  a  single  A-  or  B-family  vortex. 

The  temporal-averaged  velocity  profiles  distributed  on  a  spanwise  line  above  the  upper  surface 
of  the  delta  wing  on  a  cross  section  at  x  =  0.87L  is  shown  in  Figure  4.16.  The  distance  between 
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Figure  4.13:  Power  spectrum  density  of  v!  at  location  x  =  0.87L,  y  =  0.076L,  z  —  0.0094L 


Figure  4.14:  Power  spectrum  density  of  v 9  at  location  x  =  0.87L,  y  =  0.076L,  z  =  0.0094L 

the  spanwise  line  and  the  wing  surface  is  denoted  by  the  vertical  coordinate  z  in  Figure  4.16. 
These  results  are  in  good  agreement  with  the  experiments  carried  out  by  Rieley  &  Lowson  (1998). 
As  it  was  stated  by  Rieley  &  Lowson  (1998),  the  axial  velocity  profile  indicates  the  windward 
boundary  layer  separation.  The  inflection  point  on  the  axial  velocity  profile  is  similar  to  the 
Kelvin-Helmholtz  instability  in  plane  mixing  layers.  In  Figure  4.16(b),  the  spanwise  velocity  profile 
near  the  surface  (z— 0.0002L)  changes  sign  near  y  —  0.0746L,  which  is  corresponding  to  the  re¬ 
attachment  point  of  the  secondary  vortex.  The  inflection  points  on  profiles  of  the  vertical  velocity 
component  in  Figure  4.16(c)  are  corresponding  to  the  edge  of  the  leading-edge  shear  layer.  The 
inflection  point  moves  outboard  as  the  distance  from  the  wing  surface  increases.  The  negative 
part  of  the  vertical  velocity  is  corresponding  to  the  secondary  vortex.  The  secondary  vortex  is 
still  visible  in  the  temporal-averaged  results.  The  evidence  of  vortex  shedding  has  been  removed 
by  the  temporal  average  procedure.  The  point  of  inflection  on  the  velocity  profile  is  associated 
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Figure  4.15:  Power  spectrum  density  of  w '  at  location  x  =  0.87 L,  y  —  0.076L,  z  =  0.0094L 


with  inviscid  instability.  Rayleigh’s  inflection  point  theorem,  studying  the  instability  of  a  two- 
dimensional  velocity  profile  based  on  the  linear  normal  mode  approach,  points  out  that  disturbance 
can  be  amplified  at  the  point  of  inflection.  The  physical  interpretation  of  the  theorem  was  given 
by  Lin  (1945).  The  restoring  mechanism  will  force  a  fluid  particle  displaced  vertically  in  either 
direction  to  return  to  its  starting  position.  But  at  the  point  of  inflection  the  restoring  mechanism 
is  not  present  and  disturbance  can  grow.  In  the  flow  around  the  slender  delta  wing,  the  situation  is 
more  complex.  On  the  two-dimensional  cross  section  plane,  which  is  vertical  to  the  axial  direction, 
a  strong  shear  layer  is  observed  near  the  leading-edge,  as  shown  in  Figure  4.8(a).  The  existence 
of  the  secondary  vortex  increases  the  strength  of  the  shear  layer  and  provides  more  chance  for  the 
generation  of  the  Kelvin-Helmholtz  instability. 


4.5  Conclusions 

Direct  numerical  simulation  has  been  carried  out  to  simulate  the  flow  around  a  slender  flat- plate 
delta  wing  at  12.5°  angle  of  attack.  Two  Reynolds  number  have  been  selected.  At  a  lower  Reynolds 
number  of  5  x  104,  the  flow  is  stable  and  dominated  by  a  pair  of  leading-edge  primary  vortices. 
At  a  higher  Reynolds  number  of  1.96  x  105,  small-scale  votical  structures  are  shedding  from  the 
leading-edge.  It  has  been  found  that  the  shedding  of  the  small-scale  vortical  structures  originates 
not  only  from  the  Kelvin-Helmhotz  type  instability  of  the  leading-edge  shear  layer,  but  also  from 
the  separation  of  the  secondary  vortex  form  the  wing  surface.  The  interaction  between  the  sec¬ 
ondary  vortex  and  the  Kelvin-Helmholtz  type  vortex  is  of  particular  interesting.  The  periods  of 
vortex  shedding  are  obtained  from  the  time  series  of  velocity  components.  The  distributions  of 
the  temporal-averaged  velocity  near  the  upper  surface  of  the  delta  wing  obtained  from  the  com¬ 
putational  results  agree  well  with  those  from  the  experiment  of  Rieley  &  Lowson  (1998).  But 
the  steady  small-scale  vortical  structures  observed  in  the  same  experiment  have  not  been  found  in 
the  current  simulation  results.  However,  some  unsteady  vortical  structures  were  also  observed  in 
the  experiment,  which  could  be  related  to  the  unsteady  small  vortices  found  in  the  computational 
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(a)  Profiles  of  average  axial  velocity 


(b)  Profiles  of  average  spanwise  velocity 


(c)  Profiles  of  average  vertical  velocity 


Figure  4.16:  Variations  in  the  three  components  of  temporal-averaged  velocity  at  the  leading-edge 
with  increasing  distance  from  the  wing  surface  within  on  a  cross  section  at  x  =  0.87L 

results. 
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